######### generate nice graphs of IRF, max gain and min gain
mystar <- function(coords, v=NULL, params) {
  vertex.color <- params("vertex", "color")
  if (length(vertex.color) != 1 && !is.null(v)) {
    vertex.color <- vertex.color[v]
  }
  vertex.size  <- 1/200 * params("vertex", "size")
  if (length(vertex.size) != 1 && !is.null(v)) {
    vertex.size <- vertex.size[v]
  }
  norays <- params("vertex", "norays")
  if (length(norays) != 1 && !is.null(v)) {
    norays <- norays[v]
  }
  
  mapply(coords[,1], coords[,2], vertex.color, vertex.size, norays,
         FUN=function(x, y, bg, size, nor) {
           symbols(x=x, y=y, bg=bg,
                   stars=matrix(c(size,size/2), nrow=1, ncol=nor*2),
                   add=TRUE, inches=FALSE)
         })
}
add.vertex.shape("star", clip=igraph.shape.noclip, plot=mystar, parameters=list(vertex.norays=5))
## merge student characteristics into melted_sum_y_diff
tmp_student_data<-student_data
attach(parameters_new)
gamma_mu_study_model<-1/gamma_mu(solution$mu_study_model)  
own_slope<-(beta_vec[2]-theta_vec[1]) - theta_vec[2]*gamma_mu_study_model

tmp_student_data$own_slope<-own_slope
tmp_student_data$shocked_person<-1:n_students
tmptmp<-merge(melted_sum_y_diff, tmp_student_data, by="shocked_person")

max_val_person<-subset(tmptmp, equilibrium=="general")[which.max(subset(tmptmp, equilibrium=="general")$value),]$shocked_person
max_val_month<-subset(tmptmp, equilibrium=="general")[which.max(subset(tmptmp, equilibrium=="general")$value),]$month

tmp_degree<-3
shock_delta<-1
own_slope_scale_1<-10/.6
own_slope_scale_2<-8
y_diff_scale<-35
star_size<-20

if (max_val_month==1) A_graph_tmp<-A_1_graph else A_graph_tmp<-A_2_graph

max_val_graph<-A_graph_tmp

dist_from_max_val_person<-sapply(get.shortest.paths(A_graph_tmp, from=max_val_person)$vpath, function(x) length(x)-1)

## max val, own slope is size
# size for females: circles
  V(max_val_graph)$size<-sqrt(4/pi*own_slope)*own_slope_scale_1 - own_slope_scale_2
# size for males: squares
  V(max_val_graph)$size[student_data$male==1]<-sqrt(own_slope[student_data$male==1])*own_slope_scale_1 - own_slope_scale_2

  V(max_val_graph)$name<-1:n_students
  V(max_val_graph)$label<-NA
  V(max_val_graph)$color<-NA
  V(max_val_graph)$color[student_data$black==1]<-"grey"
  V(max_val_graph)$color[max_val_person]<-"red"
  V(max_val_graph)$shape<-"circle"
  V(max_val_graph)$shape[student_data$male==1]<-"square"
  V(max_val_graph)$size[max_val_person]<-star_size
  V(max_val_graph)$shape[max_val_person]<-"star"
  max_val_graph_tmp<-delete.vertices(max_val_graph, dist_from_max_val_person>tmp_degree | dist_from_max_val_person<0)
  set.seed(1235)
  max_layout<-layout.kamada.kawai(max_val_graph_tmp, niter=20000, kkconst=5000, initemp=20)

  ## FIGURE 1a, left
  tikz('./output/plot_max_val_graph_own_slope.tex', standAlone=F)
    plot(max_val_graph_tmp, edge.arrow.mode=0, layout=max_layout)
  dev.off()

## max val, y_diff is size
  max_val_df<-compute_shock_comparative_statics_df(A_array_months[,,max_val_month])(solution$mu_y_model, solution$mu_study_model)(solution$s_array_months_model[,max_val_month])(shocked_person_index=max_val_person, shock_delta=shock_delta)
  max_val_vec<-pmax(max_val_df[max_val_df$iteration=="general",]$y_diff, 0)

# size for females: circles
  V(max_val_graph)$size<- sqrt(4/pi*max_val_vec)*y_diff_scale
# size for males: squares
  V(max_val_graph)$size[student_data$male==1]<- sqrt(max_val_vec[student_data$male==1])*y_diff_scale
  V(max_val_graph)$size[max_val_person]<-star_size
  V(max_val_graph)$shape[max_val_person]<-"star"
  max_val_graph_tmp<-delete.vertices(max_val_graph, dist_from_max_val_person>tmp_degree | dist_from_max_val_person<0)

  ## FIGURE 1b, left
  tikz('./output/plot_max_val_graph_y_diff.tex', standAlone=F)
    plot(max_val_graph_tmp, edge.arrow.mode=0, layout=max_layout)
  dev.off()

## min val, own slope
  min_val_person<-subset(tmptmp, equilibrium=="general")[which.min(subset(tmptmp, equilibrium=="general")$value),]$shocked_person
  min_val_month<-subset(tmptmp, equilibrium=="general")[which.min(subset(tmptmp, equilibrium=="general")$value),]$month
  if (min_val_month==1) A_graph_tmp<-A_1_graph else A_graph_tmp<-A_2_graph
  min_val_graph<-A_graph_tmp
  dist_from_min_val_person<-sapply(get.shortest.paths(A_graph_tmp, from=min_val_person)$vpath, function(x) length(x)-1)

# size for females: circles
  V(min_val_graph)$size<-sqrt(4/pi*own_slope)*own_slope_scale_1 - own_slope_scale_2
# size for males: squares
  V(min_val_graph)$size[student_data$male==1]<-sqrt(own_slope[student_data$male==1])*own_slope_scale_1 - own_slope_scale_2

  V(min_val_graph)$name<-1:n_students
  V(min_val_graph)$label<-NA
  V(min_val_graph)$color<-NA
  V(min_val_graph)$color[student_data$black==1]<-"grey"
  V(min_val_graph)$color[min_val_person]<-"red"
  V(min_val_graph)$shape<-"circle"
  V(min_val_graph)$shape[student_data$male==1]<-"square"
  V(min_val_graph)$size[min_val_person]<-star_size
  V(min_val_graph)$shape[min_val_person]<-"star"
  min_val_graph_tmp<-delete.vertices(min_val_graph, dist_from_min_val_person>tmp_degree | dist_from_min_val_person<0)
  set.seed(123)
  min_layout<-layout.fruchterman.reingold(min_val_graph_tmp, niter=100000)
  
  ## FIGURE 1a, right
  tikz('./output/plot_min_val_graph_own_slope.tex', standAlone=F)
    plot(min_val_graph_tmp, edge.arrow.mode=0, layout=min_layout)
  dev.off()

## min val, y_diff
  min_val_df<-compute_shock_comparative_statics_df(A_array_months[,,min_val_month])(solution$mu_y_model, solution$mu_study_model)(solution$s_array_months_model[,min_val_month])(shocked_person_index=min_val_person, shock_delta=shock_delta)
  min_val_vec<-pmax(min_val_df[min_val_df$iteration=="general",]$y_diff, 0)
  
# size for females: circles
  V(min_val_graph)$size<- sqrt(4/pi*min_val_vec)*y_diff_scale
# size for males: squares
  V(min_val_graph)$size[student_data$male==1]<- sqrt(min_val_vec[student_data$male==1])*y_diff_scale

  V(min_val_graph)$size[min_val_person]<-star_size
  V(min_val_graph)$shape[min_val_person]<-"star"
  min_val_graph_tmp<-delete.vertices(min_val_graph, dist_from_min_val_person>tmp_degree | dist_from_min_val_person<0)
  
  ## FIGURE 1b, right
  tikz('./output/plot_min_val_graph_y_diff.tex', standAlone=F)
  plot(min_val_graph_tmp, edge.arrow.mode=0, layout=min_layout)
  dev.off()

detach(parameters_new)
